# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix C 
# Appendix C1 Figure 4 - Influence of Racial Boundary on Logged Arrests (Standardized).
# Appendix C Table 5: Influence of Racial Boundary on Logged Misdemeanor Arrests (Standardized)
# Appendix C Table 6: Influence of Racial Boundary on Logged Felonies Arrests (Standardized)
# Appendix C2 Figure 5: Louisville: Influence of Racial Boundary on Logged Arrests by Type (Standardized).
# Appendix C Table 7: Louisville: Influence of Racial Boundary on Logged Arrests by Type (Standardized).


suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# first prepare variables by city
# read in data by city

# Load Atlanta final dataset
load("atl_final.RData")

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv

# Load Austin final dataset
load("aus_final.RData")



# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv

## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv



#### atlanta: alt race measures ####

atl_b_wht = 
  lm_robust(lmisdemeanors_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + atl_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + atl_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + atl_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_whtf = 
  lm_robust(lfelonies_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

atl_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


#### austin: alt race measures ####

aus_b_wht = 
  lm_robust(lmisdemeanors_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + aus_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + aus_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + aus_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_whtf = 
  lm_robust(lfelonies_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + aus_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + aus_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + aus_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

#### boston: alt race measures ####

bos_b_wht = 
  lm_robust(lmisdemeanors_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + bos_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + bos_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_whtf = 
  lm_robust(lfelonies_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + bos_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + bos_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + bos_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime, 
            data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

#### chicago: alt race measures ####

chi_b_wht = 
  lm_robust(lmisdemeanors_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + chi_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + chi_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + chi_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_whtf = 
  lm_robust(lfelonies_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + chi_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + chi_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + chi_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

#### milwaukee: alt race measures ####

mil_b_wht = 
  lm_robust(lmisdemeanors_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + mil_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + mil_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + mil_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_whtf = 
  lm_robust(lfelonies_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + mil_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + mil_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + mil_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


#### seattle: alt race measures ####

sea_b_wht = 
  lm_robust(lmisdemeanors_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_blk = 
  lm_robust(lmisdemeanors_sd ~ p_race_black_blv + sea_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_his = 
  lm_robust(lmisdemeanors_sd ~ p_race_hisp_blv + sea_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_asn = 
  lm_robust(lmisdemeanors_sd ~ p_race_asian_blv + sea_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_whtf = 
  lm_robust(lfelonies_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_blkf = 
  lm_robust(lfelonies_sd ~ p_race_black_blv + sea_ses_blv + lpop + p_race_black + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_hisf = 
  lm_robust(lfelonies_sd ~ p_race_hisp_blv + sea_ses_blv + lpop + p_race_hisp + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea_b_asnf = 
  lm_robust(lfelonies_sd ~ p_race_asian_blv + sea_ses_blv + lpop + p_race_asian + age_15_35_male + 
              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
            data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

#### plotting out alternative race specifications #### 

df_alt_racespec = 
  data.frame(
    
    est = c(atl_b_wht$coefficients[2],
            atl_b_blk$coefficients[2],
            atl_b_his$coefficients[2],
            atl_b_asn$coefficients[2],
            atl_b_whtf$coefficients[2],
            atl_b_blkf$coefficients[2],
            atl_b_hisf$coefficients[2],
            atl_b_asnf$coefficients[2],
            
            aus_b_wht$coefficients[2],
            aus_b_blk$coefficients[2],
            aus_b_his$coefficients[2],
            aus_b_asn$coefficients[2],
            aus_b_whtf$coefficients[2],
            aus_b_blkf$coefficients[2],
            aus_b_hisf$coefficients[2],
            aus_b_asnf$coefficients[2],
            
            bos_b_wht$coefficients[2],
            bos_b_blk$coefficients[2],
            bos_b_his$coefficients[2],
            bos_b_asn$coefficients[2],
            bos_b_whtf$coefficients[2],
            bos_b_blkf$coefficients[2],
            bos_b_hisf$coefficients[2],
            bos_b_asnf$coefficients[2],
            
            chi_b_wht$coefficients[2],
            chi_b_blk$coefficients[2],
            chi_b_his$coefficients[2],
            chi_b_asn$coefficients[2],
            chi_b_whtf$coefficients[2],
            chi_b_blkf$coefficients[2],
            chi_b_hisf$coefficients[2],
            chi_b_asnf$coefficients[2],
            
            mil_b_wht$coefficients[2],
            mil_b_blk$coefficients[2],
            mil_b_his$coefficients[2],
            mil_b_asn$coefficients[2],
            mil_b_whtf$coefficients[2],
            mil_b_blkf$coefficients[2],
            mil_b_hisf$coefficients[2],
            mil_b_asnf$coefficients[2],
            
            sea_b_wht$coefficients[2],
            sea_b_blk$coefficients[2],
            sea_b_his$coefficients[2],
            sea_b_asn$coefficients[2],
            sea_b_whtf$coefficients[2],
            sea_b_blkf$coefficients[2],
            sea_b_hisf$coefficients[2],
            sea_b_asnf$coefficients[2]),
    
    se = c(atl_b_wht$std.error[2],
           atl_b_blk$std.error[2],
           atl_b_his$std.error[2],
           atl_b_asn$std.error[2],
           atl_b_whtf$std.error[2],
           atl_b_blkf$std.error[2],
           atl_b_hisf$std.error[2],
           atl_b_asnf$std.error[2],
           
           aus_b_wht$std.error[2],
           aus_b_blk$std.error[2],
           aus_b_his$std.error[2],
           aus_b_asn$std.error[2],
           aus_b_whtf$std.error[2],
           aus_b_blkf$std.error[2],
           aus_b_hisf$std.error[2],
           aus_b_asnf$std.error[2],
           
           bos_b_wht$std.error[2],
           bos_b_blk$std.error[2],
           bos_b_his$std.error[2],
           bos_b_asn$std.error[2],
           bos_b_whtf$std.error[2],
           bos_b_blkf$std.error[2],
           bos_b_hisf$std.error[2],
           bos_b_asnf$std.error[2],
           
           chi_b_wht$std.error[2],
           chi_b_blk$std.error[2],
           chi_b_his$std.error[2],
           chi_b_asn$std.error[2],
           chi_b_whtf$std.error[2],
           chi_b_blkf$std.error[2],
           chi_b_hisf$std.error[2],
           chi_b_asnf$std.error[2],
           
           mil_b_wht$std.error[2],
           mil_b_blk$std.error[2],
           mil_b_his$std.error[2],
           mil_b_asn$std.error[2],
           mil_b_whtf$std.error[2],
           mil_b_blkf$std.error[2],
           mil_b_hisf$std.error[2],
           mil_b_asnf$std.error[2],
           
           sea_b_wht$std.error[2],
           sea_b_blk$std.error[2],
           sea_b_his$std.error[2],
           sea_b_asn$std.error[2],
           sea_b_whtf$std.error[2],
           sea_b_blkf$std.error[2],
           sea_b_hisf$std.error[2],
           sea_b_asnf$std.error[2]),
    
    pv = 
      c(atl_b_wht$p.value[2],
        atl_b_blk$p.value[2],
        atl_b_his$p.value[2],
        atl_b_asn$p.value[2],
        atl_b_whtf$p.value[2],
        atl_b_blkf$p.value[2],
        atl_b_hisf$p.value[2],
        atl_b_asnf$p.value[2],
        
        aus_b_wht$p.value[2],
        aus_b_blk$p.value[2],
        aus_b_his$p.value[2],
        aus_b_asn$p.value[2],
        aus_b_whtf$p.value[2],
        aus_b_blkf$p.value[2],
        aus_b_hisf$p.value[2],
        aus_b_asnf$p.value[2],
        
        bos_b_wht$p.value[2],
        bos_b_blk$p.value[2],
        bos_b_his$p.value[2],
        bos_b_asn$p.value[2],
        bos_b_whtf$p.value[2],
        bos_b_blkf$p.value[2],
        bos_b_hisf$p.value[2],
        bos_b_asnf$p.value[2],
        
        chi_b_wht$p.value[2],
        chi_b_blk$p.value[2],
        chi_b_his$p.value[2],
        chi_b_asn$p.value[2],
        chi_b_whtf$p.value[2],
        chi_b_blkf$p.value[2],
        chi_b_hisf$p.value[2],
        chi_b_asnf$p.value[2],
        
        mil_b_wht$p.value[2],
        mil_b_blk$p.value[2],
        mil_b_his$p.value[2],
        mil_b_asn$p.value[2],
        mil_b_whtf$p.value[2],
        mil_b_blkf$p.value[2],
        mil_b_hisf$p.value[2],
        mil_b_asnf$p.value[2],
        
        sea_b_wht$p.value[2],
        sea_b_blk$p.value[2],
        sea_b_his$p.value[2],
        sea_b_asn$p.value[2],
        sea_b_whtf$p.value[2],
        sea_b_blkf$p.value[2],
        sea_b_hisf$p.value[2],
        sea_b_asnf$p.value[2]),
    
    city = c("Atlanta",
             "Atlanta",
             "Atlanta",
             "Atlanta",
             "Atlanta",
             "Atlanta",
             "Atlanta",
             "Atlanta",
             
             "Austin",
             "Austin",
             "Austin",
             "Austin",
             "Austin",
             "Austin",
             "Austin",
             "Austin",
             
             "Boston",
             "Boston",
             "Boston",
             "Boston",
             "Boston",
             "Boston",
             "Boston",
             "Boston",
             
             "Chicago",
             "Chicago",
             "Chicago",
             "Chicago",
             "Chicago",
             "Chicago",
             "Chicago",
             "Chicago",
             
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             "Milwaukee",
             
             "Seattle",
             "Seattle",
             "Seattle",
             "Seattle",
             "Seattle",
             "Seattle",
             "Seattle",
             "Seattle"),
    
    iv = 
      paste0(rep(c('Non-White', "Black", "Hispanic", "Asian", 
                   'Non-White', "Black", "Hispanic", "Asian"), 6), " Boundary"),
    
    arrtype = rep(c(rep("Misdemeanor", 4), rep("Felony", 4)), 6)
    
    
    
  )

df_alt_racespec$city = 
  factor(df_alt_racespec$city, levels = unique(df_alt_racespec$city))
df_alt_racespec$arrtype = 
  factor(df_alt_racespec$arrtype, levels = c("Misdemeanor", "Felony"),
         labels = c("A. Misdemeanor", "B. Felony"))
df_alt_racespec$iv = 
  factor(df_alt_racespec$iv,
         levels = c("Non-White Boundary",
                    "Black Boundary", "Hispanic Boundary", "Asian Boundary"),
         labels = c("Non-White Boundary",
                    "Black Boundary", "Latino Boundary", "Asian Boundary"))

# acquiring meta-analytic estimates

wht_meta = 
  df_alt_racespec %>% 
  filter(arrtype == "A. Misdemeanor") %>% 
  filter(iv == 'Non-White Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

blk_meta = 
  df_alt_racespec %>% 
  filter(arrtype == "A. Misdemeanor") %>% 
  filter(iv == 'Black Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

his_meta = 
  df_alt_racespec %>% 
  filter(arrtype == "A. Misdemeanor") %>% 
  filter(iv == 'Latino Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

asn_meta = 
  df_alt_racespec %>% 
  filter(arrtype == "A. Misdemeanor") %>% 
  filter(iv == 'Asian Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

wht_metaf = 
  df_alt_racespec %>% 
  filter(arrtype == "B. Felony") %>% 
  filter(iv == 'Non-White Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

blk_metaf = 
  df_alt_racespec %>% 
  filter(arrtype == "B. Felony") %>% 
  filter(iv == 'Black Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

his_metaf = 
  df_alt_racespec %>% 
  filter(arrtype == "B. Felony") %>% 
  filter(iv == 'Latino Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )

asn_metaf = 
  df_alt_racespec %>% 
  filter(arrtype == "B. Felony") %>% 
  filter(iv == 'Asian Boundary') %>% 
  metagen(
    data = .,
    TE = est,
    seTE = se
  )


df_alt_racespec = 
  bind_rows(
    data.frame(
      est = c(wht_meta$TE.random, blk_meta$TE.random,
              his_meta$TE.random, asn_meta$TE.random),
      se = c(wht_meta$seTE.random, blk_meta$seTE.random, 
             his_meta$seTE.random, asn_meta$seTE.random),
      pv = c(wht_meta$pval.random, blk_meta$pval.random,
             his_meta$pval.random, asn_meta$pval.random),
      iv = c('Non-White Boundary', "Black Boundary",
             "Latino Boundary", "Asian Boundary"),
      city = 'Meta-Analysis',
      arrtype = "A. Misdemeanor"
    ),
    
    data.frame(
      est = c(wht_metaf$TE.random, blk_metaf$TE.random, 
              his_metaf$TE.random, asn_metaf$TE.random),
      se = c(wht_metaf$seTE.random, blk_metaf$seTE.random, 
             his_metaf$seTE.random, asn_metaf$seTE.random),
      pv = c(wht_metaf$pval.random, blk_metaf$pval.random, 
             his_metaf$pval.random, asn_metaf$pval.random),
      iv = c('Non-White Boundary', "Black Boundary", "Latino Boundary", "Asian Boundary"),
      city = 'Meta-Analysis',
      arrtype = "B. Felony"
    )
  ) %>% 
  bind_rows(., df_alt_racespec)


df_alt_racespec$city = 
  factor(df_alt_racespec$city, levels = unique(df_alt_racespec$city),
         labels = c('Meta-\nAnalysis',
                    unique(df_alt_racespec$city)[2:length(unique(df_alt_racespec$city))]))

df_alt_racespec$arrtype = 
  factor(df_alt_racespec$arrtype, levels = c("A. Misdemeanor", "B. Felony"),
         labels = c("A. Misdemeanor", "B. Felony"))

df_alt_racespec$iv = 
  factor(df_alt_racespec$iv,
         levels = c("Non-White Boundary", "Black Boundary",
                    "Latino Boundary", "Asian Boundary"),
         labels = c("Non-White Boundary", "Black Boundary",
                    "Latino Boundary", "Asian Boundary"))

df_alt_racespec %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est,
                 col = iv),
             size = 2,
             position = position_dodge(.5)) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = iv),
                size = .4,
                position = position_dodge(.5),
                width = 0) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = iv),
                size = .8,
                position = position_dodge(.5),
                width = 0) + 
  facet_wrap(~arrtype) + 
  labs(x = "City", 
       y = "Boundary Coefficient",
       col = "Boundary Variable") + 
  scale_color_grey(start = 0, end = .6) + 
  theme_tufte(base_size = 10) + 
  theme(axis.text.x = element_text(size = 6))

ggsave(plot = last_plot(), width = 8, height = 2.5, 
       filename = 'metabrace.jpeg')



### Misdemeanor Arrests ###
# run regressions for each city 

# atlanta
atl.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                             lviolentcrime + phys_bound + com_density, 
                           data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
# austin
aus.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                             lviolentcrime + phys_bound + com_density, 
                           data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
# boston
bos.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                             phys_bound + com_density, 
                           data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# chicago
chi.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                             lviolentcrime + phys_bound + com_density, 
                           data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# milwaukee
mil.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                             lviolentcrime + phys_bound + com_density, 
                           data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# seattle
sea.model.misd = lm_robust(lmisdemeanors_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                             hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                             lviolentcrime + phys_bound + com_density, 
                           data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


# create latex table
texreg(l = list(atl.model.misd, aus.model.misd, bos.model.misd, chi.model.misd, mil.model.misd, sea.model.misd),
       file = "regtable.mainres1.misd.tex",
       include.ci = FALSE,
       custom.coef.map = list("p_race_white_blv" = "Boundary (White)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Racial Boundary on Logged Misdemeanor Arrests (Standardized)",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Misdemeanor Arrests" = 1:6),
       label = "table:misd_mainfx",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


### Felony Arrests ###
# run regressions for each city 

# atlanta
atl.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                            lviolentcrime + + phys_bound + com_density, 
                          data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
# austin
aus.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                            lviolentcrime + phys_bound + com_density, 
                          data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
# boston
bos.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                            phys_bound + com_density, 
                          data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# chicago
chi.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                            lviolentcrime + phys_bound + com_density, 
                          data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# milwaukee
mil.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                            lviolentcrime + phys_bound + com_density, 
                          data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# seattle
sea.model.fel = lm_robust(lfelonies_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                            hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                            lviolentcrime + phys_bound + com_density, 
                          data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


# create latex table
texreg(l = list(atl.model.fel, aus.model.fel, bos.model.fel, chi.model.fel, mil.model.fel, sea.model.fel),
       file = "regtable.mainres1.fel.tex",
       include.ci = FALSE,
       custom.coef.map = list("p_race_white_blv" = "Boundary (White)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Racial Boundary on Logged Felonies Arrests (Standardized)",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Felony Arrests" = 1:6),
       label = "table:fel_mainfx",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))

### Lousiville by Arrest Type ###
# create plot with robustness values for louisville for against persons, against property, and against society
# run regressions for louisville for against persons, against property, and against society



# against persons
lou.model.persons = lm_robust(lperson_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, 
                              data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

#  against property arrests
lou.model.property = lm_robust(lproperty_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                 lviolentcrime + phys_bound + com_density, 
                               data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)
# against society
lou.model.society = lm_robust(lsociety_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, 
                              data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# now to plot the coefficients
# create df of lousiville coefficients for each model
# starting with persons arrests
coef.lou.1 <- tidy(lou.model.persons)

fit_cis_95 <- confint(lou.model.persons, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(lou.model.persons, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.lou.1 <- bind_cols(coef.lou.1, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef for against persons arrests
lou_boundary_persons <- coef.lou.1[1,]


## now for against property arrests
coef.lou.2 <- tidy(lou.model.property)

fit_cis_95 <- confint(lou.model.property, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(lou.model.property, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.lou.2 <- bind_cols(coef.lou.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef for violent arrests
lou_boundary_property <- coef.lou.2[1,]

## finally against society arrests
coef.lou.3 <- tidy(lou.model.society)

fit_cis_95 <- confint(lou.model.society, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(lou.model.society, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.lou.3 <- bind_cols(coef.lou.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef for nonviolent arrests
lou_boundary_society <- coef.lou.3[1,]

# okay so let's put them all together
## combine white boundary coefs for each of the DVS
lou_coeffs <- rbind(lou_boundary_persons, lou_boundary_property, lou_boundary_society)

# for sens analysis
# against persons
lou.model.persons.sens = lm(lperson_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                              lviolentcrime + phys_bound + com_density, 
                            data = lou_fin, subset = pop > 0)
lou.model.property.sens = lm(lproperty_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                               hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, 
                             data = lou_fin, subset = pop > 0)
lou.model.society.sens = lm(lsociety_sd ~ p_race_white_blv + ses_blv + lpop + p_race_white + age_15_35_male + 
                              hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                              lviolentcrime + phys_bound + com_density, 
                            data = lou_fin, subset = pop > 0)

# robustness value exercise 
library(sensemakr)
sense_out1 = 
  sensemakr(lou.model.persons.sens,
            treatment = "p_race_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)
summary(sense_out1)

sense_out2 = 
  sensemakr(lou.model.property.sens,
            treatment = "p_race_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

summary(sense_out2)

sense_out3 = 
  sensemakr(lou.model.society.sens,
            treatment = "p_race_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

summary(sense_out3)

rvdf = data.frame(
  
  rv_val = c(as.numeric(sense_out1$sensitivity_stats$rv_q),
             as.numeric(sense_out2$sensitivity_stats$rv_q),
             as.numeric(sense_out3$sensitivity_stats$rv_q)),
  
  variable = lou_coeffs$Variable,
  
  bound = c("451x Violent Crime", "3030x Violent Crime", "2600x Violent Crime")
  
)



lou_coeffs$se = (lou_coeffs$conf.high_95 - lou_coeffs$Coefficient) / 1.96

library(meta)

out_meta = metagen(
  TE = lou_coeffs$Coefficient,
  seTE = lou_coeffs$se
)

df_meta = data.frame(
  Variable = "Meta-analysis",
  Coefficient = out_meta$TE.random,
  df = out_meta$df.Q,
  outcome = "total_larrests_sd",
  se = out_meta$seTE.random,
  conf.low_95 = out_meta$TE.random - 1.96 * out_meta$seTE.random,
  conf.high_95 = out_meta$TE.random + 1.96 * out_meta$seTE.random,
  conf.low_90 = out_meta$TE.random - 1.645 * out_meta$seTE.random,
  conf.high_90 = out_meta$TE.random + 1.645 * out_meta$seTE.random
)

lou_coeffs = bind_rows(lou_coeffs, df_meta)
# lou_coeffs$Variable = 
factor(lou_coeffs$Variable, levels = lou_coeffs$Variable)


## make xaxis labels
Plot1labels <- c("Against Persons","Against Property","Against Society")

lou_arrest_types_plot <- 
  ggplot(lou_coeffs, 
         aes(x = outcome, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 1) +
  geom_point(aes(x = outcome, 
                 y = Coefficient)) + 
  geom_errorbar(aes(x = outcome,
                    ymin = conf.low_90,
                    ymax = conf.high_90),
                width = 0,
                size = .8) +
  geom_errorbar(aes(x = outcome,
                    ymin = conf.low_95,
                    ymax = conf.high_95),
                width = 0,
                size = .4) +
  annotate("text",
           label =  paste0("RV: ", round(rvdf$rv_val[1:3], 2),
                           "\n", rvdf$bound[1:3]),
           y = lou_coeffs$Coefficient[1:3],
           x = c(1.45, 2.45, 3.45),
           size = 4,
           family = "Times") +
  geom_hline(yintercept = lou_coeffs$Coefficient[4], 
             colour = gray(1/2), lty = 2,
             alpha = .4,
             size = .2) + 
  xlab("Type of Arrests") + ylab("Racial Boundary Coefficients") + 
  scale_color_grey(start = .6, end = 0) + 
  theme_tufte(base_size = 9) + 
  scale_x_discrete(labels = c(Plot1labels, "Meta-Analysis")) + 
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 18)
  )


ggsave(width = 8, height = 6, plot = lou_arrest_types_plot, filename = "lou_arrest_types_plot.png")


# create latex table
texreg(l = list(lou.model.persons, lou.model.property, lou.model.society),
       file = "regtable.mainres2.lou.tex",
       include.ci = FALSE,
       custom.coef.map = list("p_race_white_blv" = "Boundary (White)",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "p_race_white" = "% White",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Lousiville: Influence of Racial Boundary on Logged Arrests by Type (Standardized)",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("Against Persons","Against Property","Against Society"),
       label = "table:lou_mainfx",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))



